########## re run everything
rm(list=ls())

# load packages
library(stringr)
library(ggplot2)
library(ggbreak) 
library(readxl)
library(dplyr)
library(tidyr)
library(stargazer)
library(texreg)
library(sjPlot)
library(estimatr)
library(sandwich)
library(lmtest)
library(miceadds)
library(lme4)
library(lmerTest)
library(glmnet)
library(plm)
options(scipen=999)

# load data
# set working directory to replication file
load("~Damann Kim Tavits/6. Submissions- IO/5. Replication/PostUnit.Rdata")
load("~Damann Kim Tavits/6. Submissions- IO/5. Replication/PoliticianDayUnit.Rdata")

# number of politicians
n_distinct(polday$id) # 469
n_distinct(post$id) # 469





#########################################
###### 1. H1: posts per politician ######
#########################################



####################
##### Figure 1 #####
####################

### prep data
plot_data <- polday %>% group_by(woman, date) %>%
  summarise(date = date,
            n = sum(num.post),
            n.pol = n_distinct(id[num.post!=0]),
            ppp = n/n.pol,
            days = date - as.Date("2022-02-24")) %>% distinct()

f1 <- plot_data %>% 
  ggplot(aes(x = days, y = ppp, color=woman, shape=woman, linetype = woman)) +
  geom_point(size=2.5) + xlim(-100, 100) +
  xlab("Days Since Invasion") + ylab("Posts per Politician") + 
  geom_vline(xintercept =0, linetype="dashed", color="red") +
  geom_smooth(data=subset(plot_data, days < 0), method="lm") + 
  geom_smooth(data=subset(plot_data, days > 0), method="lm") +
  theme_classic() +theme(text = element_text(size = 15)) + 
  scale_linetype_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(3, 2)) +
  scale_color_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c("#E69F00", "#56B4E9")) +
  scale_shape_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(16, 17)) 

f1

pdf("f1.pdf")
print(f1)
dev.off()



###############
### Table 1 ###
###############

# create a pdataframe for panel regressions
panel_plm <- pdata.frame(polday, index = c("id", "date"))
mod1plm <- plm(num.post ~ daysince1 + invasion*woman + daysince2, data=panel_plm,
               model="random", index = "id")
nwse <- coeftest(mod1plm, vcov.=vcovNW)[,2]
pval <- coeftest(mod1plm, vcov.=vcovNW)[,4]


screenreg(mod1plm, override.se = nwse, override.pvalues = pval,
       custom.model.names = "Posts per politician",
       custom.coef.map = list("invasion1"="Invasion",
                             "woman1"="Woman",
                             "invasion1:woman1" = "Invasion $\\times$ Woman",
                             "daysince1"="Time",
                             "daysince2" = "Time since invasion",
                             "(Intercept)"="Intercept"),
       #custom.gof.names = c(""),
       caption = "Effect of War on Politicians' Engagement",
       caption.above = T,
       include.rsquared = F,
       include.adjrs = F,
       include.nobs=T,
       include.rs=F,
       include.variance=F,
       include.deviance=F,
       booktabs = T,
       label = "tab:its",
       stars = c(0.1, 0.05, 0.01),
       custom.gof.names = c('$N_{Observations}$'),
       custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod1plm$ercomp[[1]][2]), digits=2),
                              "$N_{Politicians}$"=469),
       reorder.gof = c(1, 3, 2)
       )




###################################################
###### 2. H2: Role-Congruent Tone and Issues ######
###################################################

###############
### Table 2 ###
###############
# Effect of War on Positive Sentiment
pfull <- pdata.frame(post, index = c("id", "dayseq"))
mod2zet1 <- plm(zet ~ daysince1 + invasion*woman + daysince2, model="random",
                index="id",
                data=pfull)
mod2nrc1 <- plm(nrc ~ daysince1 + invasion*woman + daysince2, model="random",
                index="id",
                data=pfull)
mod2finn1 <- plm(afinn ~ daysince1 + invasion*woman + daysince2, model="random",
                 index="id",
                 data=pfull)
mod2bing1 <- plm(bing ~ daysince1 + invasion*woman + daysince2, model="random",
                 index="id",
                 data=pfull)
nwse <- list(coeftest(mod2zet1, vcov.=vcovNW)[,2],
          coeftest(mod2nrc1, vcov.=vcovNW)[,2],
          coeftest(mod2finn1, vcov.=vcovNW)[,2],
          coeftest(mod2bing1, vcov.=vcovNW)[,2])
pval <- list(coeftest(mod2zet1, vcov.=vcovNW)[,4],
          coeftest(mod2nrc1, vcov.=vcovNW)[,4],
          coeftest(mod2finn1, vcov.=vcovNW)[,4],
          coeftest(mod2bing1, vcov.=vcovNW)[,4])

screenreg(list(mod2zet1, mod2nrc1, mod2finn1, mod2bing1),
          override.se = nwse,
          override.pvalues = pval,
       custom.header = list("Positive Sentiment" = 1:4),
       custom.model.names = c("SYUZHET", "NRC", "AFINN", "BING"),
       custom.coef.map = list("invasion1" = "Invasion",
                              "woman1" = "Woman",
                              "invasion1:woman1" = "Invasion $\\times$ Woman",
                              "daysince1" = "Time",
                              "daysince2" = "Time since invasion",
                              "(Intercept)" = "Intercept"),
       caption = "Effect of War on Positive Sentiment",
       caption.above = T,
       include.rsquared = F,
       include.adjrs = F,
       include.nobs=T,
       include.rs=F,
       include.variance=F,
       include.deviance=F,
       label = "tab:sentiment",
       booktabs = T,
       stars=c(0.1, 0.05, 0.01),
       custom.gof.names = c('$N_{Observations}$'),
       custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=c(round(sqrt(mod2zet1$ercomp[[1]][2]), digits=2),
                                                              round(sqrt(mod2nrc1$ercomp[[1]][2]), digits=2),
                                                              round(sqrt(mod2finn1$ercomp[[1]][2]), digits=2),
                                                              round(sqrt(mod2bing1$ercomp[[1]][2]), digits=2)),
                              "$N_{Politicians}$"=c(469, 469, 469, 469)),
       reorder.gof = c(1, 3, 2))


####################
##### Figure 2 #####
####################

### prep data
plot_data <- post %>% ungroup() %>% select(woman, label, invasion) %>% 
  group_by(woman, invasion) %>%
  mutate(n.total = n()) %>% group_by(woman, invasion, label) %>%
  mutate(n = n(),
         prop = n/n.total,
         se = sqrt(prop*(1-prop)/n.total),
         ci.min = prop - 1.96*se,
         ci.max = prop + 1.96*se) %>% unique() %>% 
  mutate(across(c(prop, se), round, 5)) %>% group_by(woman, label) %>%
  mutate(diff = prop[invasion==1] - prop[invasion==0],
         diff.se = sqrt(se[invasion==1]^2 + se[invasion==0]^2))

### plot changes in proportions for social aid, morale boost, security
f2 <- plot_data %>% filter(label %in% c("aid", "pride", "sec")) %>%
  ggplot(aes(x=label, y=diff, fill=woman)) + geom_bar(stat = "identity", position="dodge") +
  labs(x="Topic Labels", y="Change in Proportions", fill="Politician's \n Gender") + 
  scale_x_discrete(labels = c("Social aid", "Morale boost", "Security")) + 
  #ylim(min = -.1, max = .1) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9"), labels=c("Men", "Women")) +
  geom_hline(yintercept=0, linetype=2, color="black") + 
  geom_errorbar(aes(ymin=diff-1.96*diff.se, ymax=diff+1.96*diff.se), width=.5, position=position_dodge(1)) +
  theme_classic() +
  theme(text = element_text(size = 15)) +
  theme(axis.text.x= element_text(size=10, angle=45, hjust=1))
f2

pdf("f2.pdf")
print(f2)
dev.off()



###############
### Table 3 ###
###############
# Effect of War on Topic Proportions
topic_panel <- polday[complete.cases(polday), ]
topic_panel_plm <- pdata.frame(topic_panel, index = c("id", "date"))
mod21p <- plm(aid ~ daysince1 + invasion*woman + daysince2, data=topic_panel_plm,
              model="random", index=id)
mod22p <- plm(boost ~ daysince1 + invasion*woman + daysince2, data=topic_panel_plm,
              model="random", index=id)
mod23p <- plm(sec ~ daysince1 + invasion*woman + daysince2, data=topic_panel_plm,
              model="random", index=id)

nwse <- list(coeftest(mod21p, vcov.=vcovNW)[,2],
             coeftest(mod22p, vcov.=vcovNW)[,2],
             coeftest(mod23p, vcov.=vcovNW)[,2])
pval <- list(coeftest(mod21p, vcov.=vcovNW)[,4],
             coeftest(mod22p, vcov.=vcovNW)[,4],
             coeftest(mod23p, vcov.=vcovNW)[,4])

screenreg(list(mod21p, mod22p, mod23p),
          override.se = nwse, override.pvalues = pval,
       custom.header = list("Topic proportion"=1:3),
       custom.model.names = c("Social aid",
                              "Morale boost",
                              "Security"),
       custom.coef.map = list("invasion1" = "Invasion",
                              "woman1" = "Woman",
                              "invasion1:woman1" = "Invasion $\\times$ Woman",
                              "daysince1" = "Time",
                              "daysince2" = "Time since invasion",
                              "(Intercept)" = "Intercept"),
       caption = "Effect of War on Topic Proportions",
       caption.above = T,
       include.rsquared = F,
       include.adjrs = F,
       include.nobs=T,
       include.rs=F,
       include.variance=F,
       include.deviance=F,
       label = "tab:3topic",
       booktabs = T,
       stars=c(0.1, 0.05, 0.01),
       custom.gof.names = c('$N_{Observations}$'),
       custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=c(round(sqrt(mod21p$ercomp[[1]][2]), digits=2),
                                                              round(sqrt(mod22p$ercomp[[1]][2]), digits=2),
                                                              round(sqrt(mod23p$ercomp[[1]][2]), digits=2)),
                              "$N_{Politicians}$"=c(469, 469, 469)),
       reorder.gof = c(1, 3, 2))


#########################################
############# 3. Reactions ##############
#########################################


####################
##### Figure 3 #####
####################

# number of posts and reactions by gender
plot_data <- polday %>% drop_na() %>%
  mutate(netreactions = react_m*num.post) %>%
  group_by(woman, date) %>% 
  summarise(date = date,
            n = sum(num.post),
            sum.react = sum(netreactions),
            m.react = sum.react/n,
            days = date - as.Date("2022-02-24")) %>% distinct()


f3 <- plot_data %>% 
  ggplot(aes(x = days, y = m.react, color=woman, shape=woman, linetype = woman)) +
  geom_point(size=2.5) + xlim(-100, 100) + ylim(0, 2500) +
  xlab("Days Since Invasion") + ylab("Reactions per Post") + 
  geom_vline(xintercept =0, linetype="dashed", color="red") +
  geom_smooth(data=subset(plot_data, days < 0), method="lm") + 
  geom_smooth(data=subset(plot_data, days > 0), method="lm") +
  theme_classic() +theme(text = element_text(size = 15)) + 
  scale_linetype_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(3, 2)) +
  scale_color_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c("#E69F00", "#56B4E9")) +
  scale_shape_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(16, 17)) 

f3

pdf("f3.pdf")
print(f3)
dev.off()




###############
### Table 4 ###
###############
# Effect of War on Public Reactions
#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################

# mod3re <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
#               model="random", index="id", data = pfull)
# nwse <-coeftest(mod3re, vcov.=vcovNW)[,2]
# pval <-coeftest(mod3re, vcov.=vcovNW)[,4]

# screenreg(mod3re, override.se = nwse, override.pvalues = pval,
#        custom.model.names = "Reactions",
#       custom.coef.map =  list("invasion1"="Invasion",
#                                "woman1"="Woman",
#                                "invasion1:woman1" = "Invasion $\\times$ Woman",
#                                "daysince1"="Time",
#                                "daysince2" = "Time since invasion",
#                                "(Intercept)"="Intercept"),
#        caption = "Effect of War on Public Reactions",
#        caption.above = T,
#        include.rsquared = F,
#        include.adjrs = F,
#        include.nobs=T,
#        include.rs=F,
#        include.variance=F,
#        include.deviance=F,
#        booktabs = T,
#        label = "table:reaction",
#        stars = c(0.1, 0.05, 0.01),
#        custom.gof.names = c('$N_{Observations}$'),
#        custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3re$ercomp[[1]][2]), digits=2),
#                               "$N_{Politicians}$"=469),
#        reorder.gof = c(1, 3, 2))



